home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YMAPS.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  11KB  |  433 lines
  1. /********************************** YMAPS ************************************/
  2. /********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
  3.  
  4.  
  5.  
  6.  
  7. #include "yinclud.h"
  8.  
  9.  
  10.  
  11. /****************  MAPS and INITIALIZATIONS OF MAPS  *******************/
  12.  
  13. mapCubic() {            /* Cubic polynomial      C    Cubic in complex
  14.                    plane  rho*Z^3 + (C1+iC2)*Z^2 + (C3+iC4)*Z +
  15.                    C5+iC6 */
  16.     int     i;
  17.     double  y_temp[6];
  18.     double  X,
  19.             XX,
  20.             Y,
  21.             YY,
  22.             XXmYY,
  23.             XY,
  24.             XY2,
  25.             Fx,
  26.             Fy;
  27.     int     base;
  28.  
  29.     X = y[0];
  30.     Y = y[1];
  31.     XX = X * X;
  32.     YY = Y * Y;
  33.     XXmYY = XX - YY;
  34.     XY = X * Y;
  35.     XY2 = 2 * XY;
  36.  /* y_temp[0] = real part of polynomial:             */
  37.  /* y_temp[0] = imaginary part of polynomial:     */
  38. /*****    y_temp[0] = rho*X*(XX-3*YY) 
  39.         +C1*XXmYY -C2*XY2 
  40.         +C3*X - C4*Y
  41.         +C5;
  42.     y_temp[1] = rho*Y*(3*XX -YY)
  43.         +C2*XXmYY +C1*XY2 
  44.         +C4*X + C3*Y
  45.         +C6;            ************/
  46.  
  47.     y_temp[0] = C5;
  48.     y_temp[1] = C6;
  49.     if(C4 != 0) {
  50.         y_temp[0] -= C4 * Y;
  51.         y_temp[1] += C4 * X;
  52.     }
  53.     if(C3 != 0) {
  54.         y_temp[0] += C3 * X;
  55.         y_temp[1] += C3 * Y;
  56.     }
  57.     if(C2 != 0) {
  58.         y_temp[0] -= C2 * XY2;
  59.         y_temp[1] += C2 * XXmYY;
  60.     }
  61.     if(C1 != 0) {
  62.         y_temp[0] += C1 * XXmYY;
  63.         y_temp[1] += C1 * XY2;
  64.     }
  65.     if(rho != 0) {
  66.         y_temp[0] += rho * X * (XX - 3 * YY);
  67.         y_temp[1] += rho * Y * (3 * XX - YY);
  68.     }
  69.  
  70.  
  71.     if(XX + YY >= 10000) {    /* to prevent overflows */
  72.         y[0] = 100;
  73.         y[1] = 100;
  74.         dot = dots;
  75.         printf("trajectory is going to infinity ");
  76.         beep();
  77.     }
  78.     for(i = 0; i < num_lyap; i++) {
  79.                 /* this is used both in Lyapunov number
  80.                    computation and for the Newton method and
  81.                    orbit continuation routines */
  82.         base = lyapzero + vec_dim * i;
  83.     /* Cauchy-Riemann equations: Fx = Gy, Fy = -Gx */
  84.         Fx = rho * 3 * (XXmYY)
  85.             + 2 * (C1 * X - C2 * Y)
  86.             + C3;
  87.         Fy = -rho * 6 * XY
  88.             - 2 * (C1 * Y + C2 * X)
  89.             - C4;
  90.         y_temp[base]
  91.             = Fx * y[base]/* partial of real part w.r.t. X */
  92.             +Fy * y[base + 1];/* partial of real part wrt Y */
  93.         y_temp[base + 1]
  94.             = -Fy * y[base]/* partial of imag part wrt X */
  95.             +Fx * y[base + 1];
  96.     }
  97.     for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
  98.         y[i] = y_temp[i];
  99.     return;
  100. }
  101.  
  102. initcubic() {
  103.     vec_dim = 2;        /* the dimension of the Lyapunov vectors =
  104.                    phase space dim */
  105.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  106.                    computed <= vec_dim  */
  107.     lyapzero = 2;        /* y[lyapzero] is the zeroth coord of the
  108.                    zeroth lyapunov vector */
  109.     X_upper = 1.;        /* x scale */
  110.     X_lower = -1.2;
  111.     Y_upper = 1.1;        /* y scale */
  112.     Y_lower = -1.1;
  113.     y[eqn0 + 0] = 0.;    /* initialize x coord of y1 */
  114.     y[eqn0 + 1] = 0.;    /* initialize y coord of y1 */
  115.     setequal(0, 1, eqn1);    /* reinitializes y[] = y1 */
  116.  /* Now initialize all constants that are used in the map even if the initial
  117.     value is 0              */
  118.     rho = 1.;
  119.     C1 = 0.;
  120.     C2 = 0.;        /* even constants that are zero must be
  121.                    initialized here since otherwise they will
  122.                    be initialized = -9999. and the menu knows
  123.                    that constants so initialized should not be
  124.                    accessible from the menu */
  125.     C3 = 0.0;
  126.     C4 = 0.;
  127.     C5 =.38;
  128.     C6 = 0.;
  129.     map = mapCubic;        /* note mapCubic is defined in this file prior
  130.                    to its mention here so it does not have to
  131.                    be declared in initCubic */
  132.     return;
  133. }
  134.  
  135.  
  136. rotor() {            /*    x->x+y     then    y->Cy+rho*sin(x)   */
  137.     double  moduloAB();
  138.     double  s,
  139.             c;
  140.     int     i;
  141.     double  y_temp[6];
  142.     int     base;
  143.  
  144.     y_temp[0] = y[0] + y[1];
  145.     s = sin(y_temp[0]);
  146.     y_temp[1] = C1 * y[1] + rho * s;
  147.  
  148.     if(num_lyap > 0)
  149.         c = cos(y_temp[0]);
  150.     for(i = 0; i < num_lyap; i++) {
  151.         base = lyapzero + vec_dim * i;
  152.         y_temp[base]
  153.             = y[base]
  154.             + y[base + 1];
  155.         y_temp[base + 1]
  156.             = rho * c * (y[base] + y[base + 1])
  157.             + C1 * y[base + 1];
  158.     }
  159.     for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
  160.         y[i] = y_temp[i];
  161.  
  162.     y[0] = moduloAB(y[0], -pi, pi);
  163.     if(C1 == 1)        /* area preserving case */
  164.         y[1] = moduloAB(y[1], -pi, pi);
  165. }
  166.  
  167.  
  168. mapComplex() {            /* Tinkerbell -- this map is not analytic */
  169.     int     i;
  170.     double  y_temp[6];
  171.     double  X,
  172.             XX,
  173.             Y,
  174.             YY;
  175.     int     base;
  176.     X = y[0];
  177.     Y = y[1];
  178.     XX = X * X;
  179.     YY = Y * Y;
  180.     y_temp[0] = XX - YY + C1 * X + C2 * Y;
  181.     y_temp[1] = 2 * X * Y + C3 * X + C4 * Y;
  182. /*    if(XX + YY >= 100000000) /to prevent overflows /
  183.         {y[0] = 100;
  184.         y[1] =100;
  185.         dot = dots;
  186.         printf("trajectory is going to infinity ");
  187.         beep();
  188.         }
  189. */
  190.     for(i = 0; i < num_lyap; i++) {
  191.                 /* this is used both in Lyapunov number
  192.                    computation and for the Newton method and
  193.                    orbit continuation routines */
  194.         base = lyapzero + vec_dim * i;
  195.         y_temp[base]
  196.             = (2 * X + C1) * y[base]
  197.             + (-2 * Y + C2) * y[base + 1];
  198.         y_temp[base + 1]
  199.             = (2 * Y + C3) * y[base]
  200.             + (2 * X + C4) * y[base + 1];
  201.     }
  202.     for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
  203.         y[i] = y_temp[i];
  204.     return;
  205. }
  206.  
  207.  
  208.  
  209.  
  210.  
  211. initComplex() {
  212.     vec_dim = 2;        /* the dimension of the Lyapunov vectors =
  213.                    phase space dim */
  214.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  215.                    computed <= vec_dim  */
  216.     lyapzero = 2;        /* y[lyapzero] is the zeroth coord of the
  217.                    zeroth lyapunov vector */
  218.     X_upper =.5;        /* x scale */
  219.     X_lower = -1.3;
  220.     Y_upper = -1.6;        /* y scale */
  221.     Y_lower = 0.6;
  222.     y[eqn0 + 0] = -.72;
  223.     y[eqn0 + 1] = -.64;
  224.     setequal(0, 1, eqn1);    /* reinitializes y[] */
  225.     C1 =.9;
  226.     C2 = -.6013;
  227.     C3 = 2.0;
  228.     C4 =.5;
  229.     map = mapComplex;
  230. }
  231. initRotor() {
  232.     vec_dim = 2;        /* the dimension of the Lyapunov vectors =
  233.                    phase space dim */
  234.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  235.                    computed <= vec_dim */
  236.     lyapzero = 2;        /* y[lyapzero] is the zeroth coord of the
  237.                    zeroth lyapunov vector */
  238.  
  239.     X_upper = pi;        /* x scale */
  240.     X_lower = -pi;
  241.     Y_upper = pi;        /* y scale */
  242.     Y_lower = -pi;
  243.     C1 = 1;
  244.     rho = 1.;
  245.     map = rotor;        /*  31 */
  246. }
  247.  
  248. LaserMap() {            /* Ring Laser map:(Z= x+iy) Z -> rho+ C2*Z*exp{
  249.                    i[C1 - C3/(1 + |Z|**2)] } */
  250.     double  y_temp[6];
  251.     int     base,
  252.             i;
  253.     double  cs,
  254.             sn,
  255.             temp,
  256.             xn,
  257.             ynval,        /* YN is a global variable */
  258.             C3over,
  259.             real,
  260.             imag;
  261.     double  tempprime,
  262.             dfdx,
  263.             dfdy,
  264.             dgdx,
  265.             dgdy;        /* first derivatives */
  266.     xn = y[0];
  267.     ynval = y[1];
  268.     C3over = C3 / (1 + xn * xn + ynval * ynval);
  269.     temp = C1 - C3over;
  270.     cs = C2 * cos(temp);
  271.     sn = C2 * sin(temp);
  272.     real = xn * cs - ynval * sn;
  273.     imag = xn * sn + ynval * cs;
  274.     y_temp[0] = rho + real;
  275.     y_temp[1] = imag;
  276.  
  277.     if(num_lyap > 0) {
  278.         tempprime = 2 * C3over * C3over / C3;
  279.         dfdx = cs - imag * tempprime * xn;
  280.         dfdy = -sn - imag * tempprime * ynval;
  281.         dgdx = sn + real * tempprime * xn;
  282.         dgdy = cs + real * tempprime * ynval;
  283.  
  284.         for(i = 0; i < num_lyap; i++) {
  285.             base = lyapzero + vec_dim * i;
  286.             y_temp[base]
  287.                 = dfdx * y[base]
  288.                 + dfdy * y[base + 1];
  289.             y_temp[base + 1]
  290.                 = dgdx * y[base]
  291.                 + dgdy * y[base + 1];
  292.         }
  293.     }
  294.  
  295.     for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
  296.         y[i] = y_temp[i];
  297. }
  298.  
  299. LaserInverseMap() {        /* Ring Laser map:(Z= x+iy) Z -> rho+ C2*Z*exp{
  300.                    i[C1 - C3/(1 + |Z|**2)] } */
  301.     double  C2over,
  302.             psi,
  303.             ypsi,
  304.             xpsi,
  305.             dxpsi,
  306.             dypsi;
  307.     double  arg,
  308.             dxarg,
  309.             dyarg;
  310.     double  y_temp[6];
  311.     int     base,
  312.             i;
  313.     double  cs,
  314.             sn,
  315.             xn,
  316.             ynval;
  317.     double  dfdx,
  318.             dfdy,
  319.             dgdx,
  320.             dgdy;        /* first derivatives */
  321.  
  322.  
  323.     if(C2 == 0) {
  324.         PRINT "IMPOSSIBLE VALUE: C2 = 0; no inverse exists\n");
  325.         return;
  326.     }
  327.     xn = y[0];
  328.     ynval = y[1];
  329.     C2over = 1 / C2;
  330.     xpsi = (xn - rho) * C2over;
  331.     ypsi = ynval * C2over;
  332.     psi = 1 + xpsi * xpsi + ypsi * ypsi;
  333.     arg = C1 - C3 / psi;
  334.     cs = C2over * cos(arg);
  335.     sn = C2over * sin(arg);
  336.     y_temp[0] = (xn - rho) * cs - ynval * sn;
  337.     y_temp[1] = ynval * cs + (xn - rho) * sn;
  338.  
  339.  
  340.  
  341.     if(num_lyap > 0) {
  342.         dxpsi = 2 * (xn - rho) * C2over;
  343.         dypsi = 2 * ynval * C2over;
  344.         dxarg = -C3 * dxpsi / (psi * psi);
  345.         dyarg = -C3 * dypsi / (psi * psi);
  346.         dfdx = cs - ((xn - rho) * sn + ynval * cs) * dxarg;
  347.         dfdy = ((-xn + rho) * sn - ynval * cs) * dyarg - cs;
  348.         dgdx = (-ynval * sn + (xn - rho) * cs) * dxarg - sn;
  349.         dgdy = cs + (-ynval * sn + (xn - rho) * cs) * dyarg;
  350.  
  351.         for(i = 0; i < num_lyap; i++) {
  352.             base = lyapzero + vec_dim * i;
  353.             y_temp[base]
  354.                 = dfdx * y[base]
  355.                 + dfdy * y[base + 1];
  356.             y_temp[base + 1]
  357.                 = dgdx * y[base]
  358.                 + dgdy * y[base + 1];
  359.         }
  360.     }
  361.  
  362.     for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
  363.         y[i] = y_temp[i];
  364. }
  365.  
  366. initLaser() {
  367.     vec_dim = 2;        /* the dimension of the Lyapunov vectors =
  368.                    phase space dim */
  369.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  370.                    computed <= vec_dim */
  371.     lyapzero = 2;        /* y[lyapzero] is the zeroth coord of the
  372.                    zeroth lyapunov vector */
  373.  
  374.     X_upper = 3;        /* x scale */
  375.     X_lower = -1;
  376.     Y_upper = 4.5;        /* y scale */
  377.     Y_lower = -4.0;
  378.     C2 = 0.9;
  379.     rho = 1.;
  380.     C3 = 6.;
  381.     C1 =.4;
  382.     map = LaserMap;        /*  50  */
  383. }
  384.  
  385. bob() {            /* this routine defines the map */
  386.     int     rand();    /* a Microsoft random integer generator */
  387.     double  temp0,
  388.             temp1,
  389.             c,
  390.             s,
  391.             angle;
  392.  
  393.     if(rand () < C2 * 32768.) {
  394.                 /* rand() returns a random integer between 0
  395.                    and 32768, so C2 should be a number between
  396.                    0 and 1, so that the first process will be
  397.                    carried out C2 of the time */
  398.         y[0] = 2.- y[0];
  399.         y[1] = -y[1];
  400.     }
  401.     else {
  402.         angle = C1 * twopi / 360.;
  403.         c = cos(angle);
  404.         s = sin(angle);
  405.         temp0 = rho * (c * y[0] + s * y[1]);
  406.         temp1 = rho * (-s * y[0] + c * y[1]);
  407.         y[0] = temp0;
  408.         y[1] = temp1;
  409.     }
  410. }
  411.  
  412. init_bob() {            /* this process initializes the constants and
  413.                    with the first line tells the program the
  414.                    name of the routine that defines the map;
  415.                    since C1, C2, and rho are used in bob(),
  416.                    htey must be initiallized in this routine;
  417.                    if they are not, the program will not let be
  418.                    accessed interactively and they will be
  419.                    given default values, namely -9999. */
  420.     map = bob;        /* this tells the program the name of the
  421.                    routine that is evaluated; map is a pointer
  422.                    to routines and bob() must be defined in the
  423.                    file prior to this line; */
  424.     rho =.65;
  425.     C1 = 135.;
  426.     C2 =.5;
  427.  /* the following 4 determine the x and y scale */
  428.     X_lower = -3.;
  429.     X_upper = 5.;        /* x scale */
  430.     Y_lower = -3.;
  431.     Y_upper = 3.;        /* y scale */
  432. }
  433.